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Abstract:  We  calculate  the  absorption  spectrum  of  a  cluster  using  a 
computational  method  in  which  classical-like  degrees  of  freedom  for  the  nuclei 
are  described  by  Gaussian  wave  packets  while  the  valence  electrons  are  treated 
quantum  mechanically.  We  examine  the  spectral  features  in  comparison  to  an 
even  simpler  mixed  quantum-classical  model  in  which  the  nuclear  motion  is 
treated  by  purely  classical  mechanics.  Anomalous  features  (  such  as  negative 
absorption  )  in  the  absorption  spectrum  which  can  arise  from  mixed  quantum- 
classical  methods  are  examined  and  the  Gaussian  wave  packet  nuclear  dynamics 
is  found  to  substantially  reduce  these  anomalous  features.  This  method  is  applied 
to  a  two-coordinate  model  problem  in  which  exact  numerical  results  can  be 
obtained  and  we  find  that  the  method  works  fairly  well.  We  also  apply  the  method 
to  the  valence  electronic  absorption  spectrum  for  a  KXe6  cluster.  The  method  does 
not  suffer  from  the  dramatic  failure  seen  when  Xe  motion  is  treated  classically. 
The  method  is  used  to  calculate  the  vibrational  width  and  the  vibrational 
structure  of  the  electronic  absorption  spectrum. 
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L  Introduction 

There  is  an  enormous  literature  concerning  computational  methods  and 
applications  in  which  some  degrees  of  freedom  are  treated  quantum 
mechanically  and  some  are  treated  classically  or  semiclassically^1'21^  The  basic 
idea  is  that  light  particles,  such  as  electrons,  require  a  quantum  description 
while  nuclei  may  or  may  not  depending  upon  the  specific  application. 

There  is  also  a  large  literature  concerning  the  errors  made  by  treating 
some  degrees  of  freedom  classically  and  others  quantum  mechanically  113,16,22-28] 
There  is  as  yet  no  precise  criterion  that  can  be  used  to  establish  in  advance  when 
one  is  guaranteed  to  obtain  accurate  results,  but  the  existing  computational 
experience  allows  us  to  guess  the  likelihood  of  success.  The  outcome  depends  to  a 
large  extent  on  the  quantity  being  computed.  Calculations  of  the  expectation  value 
of  an  observable,  given  by  an  expression  of  the  form  ,  (  T(t)  |  O  |  T(t) )  ,  where  0  is 
an  operator  and  1  *F(t ) )  is  the  time  dependent  wave  function  of  the  system  are 
often  successful.  They  do  however  tend  to  fail  if  the  wave  function  for  the 
quantum  degrees  of  freedom  consists  of  "pieces"  localized  in  different  regions  of 
the  spaceJ11,13,16,26^  Quantities  of  the  form  <  TCt)  |  O  |  T(0)  )  which  appear  in 
calculations  of  various  rates  (  e.g.  photon  absorption  and  Raman  scattering, 
radiationless  transitions,  etc.)  are  more  difficult.  Model  calculations  have  shown, 
for  example,  that  calculations  of  the  absorption  cross  section  lead  to  unacceptable 
results,  such  as  fairly  large  negative  absorption  rates.  The  reason  for  this  failure 
was  determined®-21^  by  numerical  analysis  using  various  levels  of  approximations, 


for  a  model  for  which  exact  numerical  results  were  also  obtained.  In 
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conventional  mixed  quantum-classical  (QC)  calculations  only  the  evolution  of  the 
quantum  degrees  of  freedom  is  represented  by  a  wave  function.  The  time 
dependent  positions  of  the  particles  treated  classically  appear  parametrically  in 
the  Hamiltonian  for  the  quantum  degrees  of  freedom.  The  evolution  of  these 
classical  degrees  of  freedom  -  obtained  by  solving  an  'effective''  Newton  equation  - 
provides  a  time  dependent  force  acting  on  the  quantum  degrees  of  freedom.  For 
example,  the  quantum  particle  could  be  an  electron  and  the  classical  particles, 
nuclei.  If  the  nuclei  oscillate  at  a  certain  frequency,  the  force  they  exert  on  the 
electron  will  oscillate  at  this  frequency,  and  one  may  expect  that  this  oscillating 
force  will  correctly  generate  the  vibrational  structure  in  the  spectrum.  The 
results^2 ^  do  not  bear  out  this  expectation.  The  reason  is  simple.  In  the  quantum 
theory,  the  spectral  information  is  contained  in  the  time  evolution  of  the  overlap 
integrals  between  the  initial  nuclear  wave  function  and  the  same  wave  function  at 
time  t.  The  QC  theory  does  not  provide  a  wave  function  for  the  nuclear  wave 
functions  and  therefore  this  overlap  is  set  equal  to  one.  This  error  is  compounded 
by  the  fact  that  the  nuclei  act  on  the  electron  with  a  time  dependent  force.  If  the 
nuclei  were  not  allowed  to  move,  setting  the  nuclear  overlap  to  one  will  merely  kill 
the  vibrational  bands  in  the  spectrum;  moving  them  classically  brings  about 
negative  rates. 

Metiu  and  Haug  ^  ,  who  pointed  out  this  unpleasant  feature,  suggested 
a  model  that  will  remedy  the  situation  without  adding  much  to  the  computational 
burden:  to  treat  the  nuclear  motion  by  Heller's  Gaussian  Wave  Packet  (GWP) 
method. f29a]  This  is  known  to  give  good  results  in  short  time  calculations  and  to 
fail  for  long  times.  Procedures  that  will  give  better  long  time  results^291^  tend  to 
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become  computationally  expensive,  defeating  the  main  reason  for  the  GWP 
approximation.  This  impasse  is  resolved  by  an  important  observation  made  by 
Nitzan  (private  communication),  in  the  context  of  the  theory  of  radiationless 
transitions,  which  is  also  applicable  to  absorption  rate  calculations  considered 
here.  In  our  case  we  must  calculate  the  evolution  of  the  system  in  the  "promoted 
state"  *-30^  ,  r  |  ¥(0)  > ,  where  r  is  the  electron  position  and  |  T(0) )  is  the  ground 
state  of  the  system.  The  nuclei,  which  started  out  at  their  equilibrium  positions, 
corresponding  to  the  state  |  T(0) ) ,  are  out  of  equilibrium  with  respect  to  the  state  r 
I  VF(0)  )  .  As  a  result,  the  nuclei  will  start  to  move  as  soon  as  the  electron  is 
"promoted"  and  the  overlap  integral  of  their  current  wave  function  with  their 
initial  wave  function  will  start  going  down.  If  the  electron  interacts  with  N 
nuclear  degrees  of  freedom,  we  will  have  a  product  of  N  overlap  integrals,  all 
becoming  less  than  one.  While  each  overlap  may  change  little  (  i.e.  from  1.0  to  0.8 
)  the  product  will  plummet  rapidly.  Thus  for  systems  with  many  degrees  of 
freedom  we  only  need  a  short  time  calculation  to  determine  the  vibrational 
broadening  of  the  spectrum.  During  this  short  time,  the  nuclei  cannot  sample 
much  of  the  potential  energy  surface,  which  can  then  be  represented  by  a  local 
harmonic  approximation.  .  Because  of  this  a  GWP  calculation  of  the  nuclear 
overlap  should  be  very  accurate. 

In  this  article  we  explore  numerically  how  this  works  for  two  cases.  In 
both  cases,  there  is  a  single  valence  electron  treated  quantum  mechanically  while 
the  nuclei  are  treated  by  GWP  -  and  we  denote  this  combined  method  QGWP.  The 
first  case  is  a  two-coordinate  test  model  mimicking  an  alkali-helium  dimer  that 
we  have  previously  used^2^  for  the  examination  of  other  QC  methods.  In  this 
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case,  because  of  the  simplicity  of  the  model,  we  can  examine  exact  quantum 
mechanical  and  also  exact  quantum  Hartree  calculations  as  benchmarks  for 
comparison  with  the  QGWP  method.  The  second  case  we  examine  is  an 
application  to  a  3-dimensional  cluster  with  realistic  potentials  mimicking  KXeg. 
In  this  cluster,  the  three  dimensional  motion  of  the  valence  electron  is  treated 
quantum  mechanically  while  the  nuclear  motion  is  treated  by  the  GWP 
approximation.  The  initial  condition  is  the  quantum  mechanical  ground  state  for 
the  electron  and  the  GWP  approximation  for  the  ground  state  of  the  nuclei.  To 
calculate  the  absorption  cross  section,  the  electron  is  promoted  and  propagated  by 
a  Fast  Fourier  Transform  (FFT)  method^33^  while  the  nuclear  wave  function 
evolves  according  to  Heller's  GWP  method. 

In  Section  II,  we  outline  the  theory  behind  these  calculations.  In  Section 
III  the  two-coordinate  test  model  is  presented  and  in  Section  IV  the  KXe6  cluster 
is  examined. 


IL  THEORY 


ELA.  Introduction 

We  consider  a  cluster  containing  a  single  valence  electron  (located  at  r), 
the  parent  ion  core  with  charge  +1  and  N  "solvent”  atoms  (  located  at  R  =  { Ra  , 
a=l,2,3...,N+l  }  ).  The  corresponding  momenta  are  p  and  P  respectively.  The 
Hamiltonian  is 

H  =  Ha(R  ,P)  +  He(r,p;R).  (II.  1) 

Ha  ,  describing  the  atomic  degrees  of  freedom  (i.e.  the  motion  of  the  ion 
core  and  the  solvent  atoms),  is  given  by 

N+l 

Ha  =  X  Pa2  /  (2Ma)  +  Vg(R) ,  (II. 2) 

a=l 

where  Va(R)  is  the  sum  of  the  ion-solvent  and  the  solvent-solvent  interaction 
energies. 

He  ,  describing  the  valence  electron  degrees  of  freedom,  is  given  by 

He=  p2/(2m)  +  Vea(r3)  (H-3) 

containing  the  kinetic  energy  of  the  electron  plus  the  interaction  energy,  Vea  , 
which  is  a  sum  of  the  electron-ion  and  the  electron-solvent  interaction  energies. 
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IT.B.  Absorption  Spectrum 

To  calculate  the  absorption  cross  section  a  (to)  for  the  cluster  configuration 
T  ,  we  use  a  slight  variation  of  Heller’s  formula^ 


o(co)  ~  co Kef  dt  exp(icot)  C(t) .  (II. 4) 

Jo 

The  overlap  integral  C(t)  is  given  by 

C(t)=  exp(iEgt/h)  ( Tp  (0)  lTp(t)),  (II. 5) 

where  Eg  is  the  total  ground  state  energy  of  the  system  and  where  the  promoted 
state  at  time  zero,  I  Tp  (0)  > ,  is 

I  H'p  (0)  >  =  -e  E*r  i  Tg(0))  (II.6) 

In  Eq.(II.6)  1  'Fg  (0) )  is  the  electronic  and  nuclear  ground  state,  -er  is  the  dipole 
operator  and  e  is  the  electric  field  polarization  direction.  The  promoted  state  at 
time  t  is 


I  'Fp(t) )  =  exp(-iHt/h  )  I  Tp  (0) ) .  (II. 7) 

This  formulation  provides  the  absorption  spectrum  to  all  the  excited  states 
(including  the  photo-emission  spectrum)  in  one  calculation.  The  computations 
needed  to  evaluate  Eq.(II.5)  are  made  possible  by  using  pseudopotentials  to 
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describe  the  electron-ion  and  electron-solvent  interactions  resulting  in  an  one 
electron  problem,  and  by  approximating  the  nuclear  dynamics  in  some  manner. 

ILC.  A  summary  of  the  time  dependent  Hartree  approximation 

In  all  QC  methods  the  division  of  the  degrees  of  freedom  into  classical  and 
quantum  ones  is  made  by  using  a  Hartree  product  wave  function  approximation. 
For  calculating  the  spectrum,  we  use  a  Hartree  product  for  the  promoted  state 
wave  function  [see  Eq.(II.5)  ] 

<  r,R  1  Tp  )  =  (rl-erl  bg>  ( Rl  Xg)exp(irj/h)  =  ( r  I  $p)  (  Rl  xg>exp(in/h).  (II. 8) 

Since  our  interest  here  is  in  electric  dipole  spectroscopy,  the  dipole  operator  is  -er, 
and  therefore  only  the  ground  state  electronic  wave  function,  1 4>  )  ,  is  'promoted'' 
to,  1  <j>p)  ,  while  the  ground  state  nuclear  wavefunction,  I  xg>  ,  is  unaffected.  In  Eq. 
(II. 8),  as  is  typical  in  time-dependent  Hartree  (tdH)  approximations,  a  time 
dependent  phase  term,  exp(iq(t)/h)  ,  is  separated  out  to  simplify  the  form  of  the 
tdH  equations  as^32^ 

ih|  I  ^p)  =  [  p2/2m  +  (Xg  I  Vea(r,R)  1  %g)  /  (xg  I  Xg)  3  1 0p>,  UI.9) 

ih|-  lxg>  =  [p2/2M  +  Va(R)  +  {»p|Vea(r,R)|1t,p)/(0pl9p)!  I*.).  (11.10) 

and 

I  ’Ig  -  1  Vea<r.R> 1  XgV  '  <Xg  I  Xg>  <*p  'V  1  n- 1 1  > 


The  Hartree  product  wave  function  for  the  promoted  state  applied  to  the 
correlation  function  of  Eq.(II.5)  leads  to  a  “Hartree1'  correlation  function 

Ch(t)  =  exp(iEgt/h)  (<t>p(0)  I  <f>p(t))  (xg(0)  I  xg(t)>  exp(i[r|g(t)-ng(0)]/h) , 

(11.12) 

with  a  Hartree  ground  state  energy  given  by 

Eg  =  <4>e  I  P2/2m  1 4,g>  /  Wg  1 4g>  +  <xg  I  P2/2M+Va(R)  I  Xg>  /  <Xg  I  Xg>  + 

<<?gXg  I  Vea(r,R)  I  Xg*6>  /  (Xg  I  Xg)  <  W  •  (II.  1 3) 

We  note  that  due  to  the  broken  correlation  in  the  Hartree  product  state, 
that  the  effective  potentials  (the  averages  over  Vga  )  ir  Eq.(II.9-ll)  become  time 
dependent  send  can  give  rise  to  time  varying  energy  levels  in  the  system.  This 
variation  of  the  energy  levels  in  time  can  result  in  serious  problems  in  the 
absorption  spectrum:  unphysical  negative  amplitudes  can  appear  as  examined  in 
Ref.  21  .  Therefore  any  method  that  begins  with  the  Hartree  separation,  such  as 
the  one  we  use  here,  will  be  susceptible  to  this  problem.  A  failure  to  include  the 
amplitude  and  phase  information  arising  from  the  classical-like  degrees  of 
freedom  will  tend  to  magnify,  possibly  to  a  catastrophic  extent  (  Ref.  21  ),  the 
spectroscopic  anomalies  that  can  arise  even  at  the  Hartree  level.  Our  quest  is  to 
minimize  these  anomalous  effects  by  building  an  approximate,  easy  to  evaluate 
wave  function.  By  using  the  GWP  method  we  can  generate  the  quantum 
amplitude  and  phase  information  associated  with  the  classical-like  nuclear 
degrees  of  freedom,  which  is  needed  for  the  absorption  cross  section  calculation. 
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ILD  The  QGWP  method 

The  wave  function  ,  xg(R,t)  ,  representing  the  "classical-like"  nuclear 
degrees  of  freedom,  R  ,  and  the  integrals  involving  xg  will  be  approximated 
because  of  the  computational  complications  arising  from  the  many  dimensional 
character  of  xg-  If  we  examine  the  Hartree  equation  (II.  10)  for  %g(R,t)  we  see  that 
the  effective  Hamiltonian  for  these  "classical-liKe"  degrees  of  freedom  is 

Heff  =  P2/2M+Veff  (11.14) 

Veff=  Va(R)  +  <(j>p  1  Vea(r,R) !  <fjp)  /  (bp  1  $p) .  (11.15) 

Moving  to  the  classical  limit  we  take  Eq.(II.14)  as  the  effective  classical 
Hamiltonian  for  the  atoms  in  the  cluster.  Classical  trajectories  for  the  effective 
classical  Hamiltonian  in  terms  of  the  Cartesian  coordinates  and  momenta  are 


given  by 

A  Pa  =  '  VeMR)  ■  <n,16) 

dt  dKa 

d_Ra  =  Pa/Ma.  (II.17) 

dt 

In  the  QC  dynamic  mixture  we  also  make  approximations  for  the 
integrals  occurring  in  Eq.(II.9,  11  and  13)  over  the  R  coordinates  of  the  forms: 

<Xg  I  Vea(r,R)  I  xg)  /  (xg  I  Xg)  ~  Vea(  r,  (R) ) ,  (11.18a) 

< Vp  1  Vea^R> 1  Xgbp>  /  (Xg  1  Xg)  (bp  I V  ~  <*p  I  Vea(  r,  (R) )  I  <j>p)  /  <bp  I  0P) ,  (II.18b) 

(bgXg  1  Vea(r,R)  I  xgbg)  /  (Xg  l  Xg)  (bg  I  bg>  ~  (bg  I  Vea(  r,  (R> )  I  bg)  /  (bg  I  Og) ,  (II.  18c) 
(Xg  I  P2/2M+Va(R)  I  xg)  /  (Xg  I  Xg)  ~  (P2)/2M  +  Va(  <R) ) ,  (II.  18d) 

where  in  the  classical  limit 

OR)  =  (Xg  I R I  Xg>  /  <Xe I  Xg>  -  R(t) .  (II.19a) 

,P2)  =  (ZgIP2IXg)'(XgIXg)-P(t)2. 
are  represented  by  the  classical  trajectories  for  R(t)  and  P(t). 


(II. 19b) 
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The  initial  condition  of  the  duster,  for  a  specified  cluster  size,  is 
determined  by  finding  a  configuration  of  (possibly  local)  minimum  energy.  To  be 
consistent  with  the  classical  equations,  the  minimum  energy  atomic 
configuration  for  the  cluster  is  determined  with  respect  to  the  effective  classical 
potential  Veff  by  following  its  gradient  with  respect  to  a  parameter  t 

^R  =  -VVeff(R).  (11.20) 

Fluctuations  are  introduced  in  the  early  phase  of  the  minimization  algorithm  so 
as  to  avoid  the  likelihood  of  falling  into  a  local  basin  of  the  effective  potential.  The 
fluctuations  are  introduced  by  Monte  Carlo  (MC)  sampling  from  the  distribution 

exp{-(3Va(R)  }  exp{-{3  <$gl  Vea(r,R)lbg>  1 ,  (11.21) 

in  which  (3  gradually  evolves  from  a  high  temperature  (typically  1000K  for  the 
applications  below)  to  a  low  temperature  (  typically  100K),  as  various 
configurations  are  generated  by  the  MC  procedure.  Each  sampling  provides  the 
ion  core  and  the  solvent  atom  positions  Rj,  {Ra}a=1  n  corresponding  to  a  highly 
probable  configuration.  Because  the  electron  excitation  gap  is  much  larger  than 
kBT  ,  the  electron  is  in  the  ground  state  and  the  matrix  element  appearing  in  Eq. 
(11.21)  is  equal  to  the  ground  state  interaction  energy  of  the  electron  for  the  given 
nuclear  configuration.  After  this  "cooling  down"  period  we  follow  the  gradient  to 
the  minimum  effective  potential  configuration.  The  net  effect  of  choosing  the 
initial  configuration  in  this  way  is  that  if  the  quantum-classical  dynamic  scheme 
is  initiated  without  promoting  the  electron  then  nothing  happens:  the  electron 
remains  in  the  ground  state  and  the  cluster  remains  frozen  in  its  configuration, 
which  is  the  required  behavior. 

As  a  further  computational  convenience,  we  make  an  electronic  Born- 
Oppenheimer  approximation  by  ignoring  the  parametric  R  dependence  of  $  (r,t;R> 
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when  evaluating  the  gradient  in  Eq.(II.16  and  20);  i.e.  we  drop  the  term  involving 
^<J)p(r,t;R).  We  note  that  this  parametric  R  dependence  of  <j)p(r,t;R)  is  included  in 
the  MC  sampling  of  Eq(II.21)  and  if  we  established  an  initial  configuration  solely 
by  the  use  of  the  MC  sampling  function  (damped  gradually  to  0  K),  a  slightly 
different  configuration  would  result  than  that  given  by  the  gradient  method. 

From  the  local  minimum  configuration,  the  initial  condition  for  the  GWP 
approximation  for  xg(R>t)  is  determined  via  a  normal  mode  analysis  about  this 
stable  configuration.  The  function  xg(R,t)  is  then  given  by 

Xg(R,t)  =  []  gn(t)  (11.22) 

n 

where  n  runs  over  the  3N-6  vibrational  modes  of  the  N-atom  cluster.  Qn  and  con 
are  The  mass  scaled  normal  mode  coordinates,  Qn  ,  and  normal  mode 
frequencies  ,  con  ,  satisfying 

Qn  +  con2  Qn  =  0  ,  (11.23) 

are  used  to  specify  the  initial  normal  mode  wave  functions 

gn(q,t=0)  =  (con/7tft)1/4  exp[  -  q  -  Qn(0) }  2  /  2h  ]  (11.24) 

used  as  the  initial  ground  state  wave  function  for  the  n-th  mode.  Inserting 
Eq.(II.24)  into  Eq.(II.22)  gives  the  initial  nuclear  wave  function. 

The  time  evolution  of  the  gn(q,t)  wave  functions  is  performed  by  Heller's 
version  of  Gaussian  Wave  Packet  (GWP)  dynamics^'  resulting  in 

gn(q,t)  =  (co^h)174  exp{  (i/h)[  (icon/2){  q  -  Qn(t) }  2  +  nn(t)  (  q  -  Qn(t)  1  -  yn(t)  ]  j  . 


(11.25) 
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We  have  used  the  initial  conditions,  given  by  Eq.(II.24),  to  set  nn(0)  =  0  ,  yn(0)  =  0  , 
and  also  to  argue  that  the  frequencies,  con  ,  which  serve  as  time-dependent  width 
parameters  within  the  general  form  of  the  GWP  dynamics,  should  be  time 
independent  for  the  relatively  short  times  involved  in  low  resolution  spectral 
calculations  (this  is  the  frozen  Gaussian  version  of  GWP  dynamics*-2^). 

The  variable  n  is  the  mass  scaled  normal  mode  momenta  ,  i.e. 

na  =  AQa  (11.26) 

dt 

and  the  Gaussian  phase  parameter  is  determined  from 

A  Ya  =  (l/2)na2-Veff(Qa)  -  h(V2  .  (11.27) 

dt 

In  Eq.(II.27)  we  use  the  local  harmonic  approximation 

veff  <Qtt>  »  (1/2)  0Jo2  Qa2  +  vzero  (11.28) 

with 

Vzero  =  Va(  R  )  +  I  vea(  r,E  )  Hp>  /  <+p  I  ♦„)) .  (11.29) 

being  evaluated  with  <j>p(t=0)  and  R(t=0).  We  emphasize  that  because  we  need  only 
the  short  time  behavior  of  the  nuclear  overlap  integral,  the  local  harmonic 
approximation  should  work  very  well. 

We  can  interpret  the  dynamics  expressed  in  the  promoted  state 
correlation  function  by  saying  that  the  promotion  of  the  electron  from  the  ground 
to  the  promoted  state  provides  the  energy  to  set  the  cluster  into  motion.  So  both  the 
electronic  excitations  and  the  vibrational  structure  of  those  excitations  are 
included  in  the  correlation  function. 


DUE  The  trajectory  approximation 

A  first  approximation  often  used  for  condensed  phase  correlation 
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functions  involving  a  system  with  a  localized  quantum  component  surrounded  by 
more  classical-like  degrees  of  freedom  is  to  take  a  mean  field  approach  with 
respect  to  the  classical-like  coordinates.  This  approach  builds  in  the  effect  of  the 
average  perturbation  of  the  surrounding,  classical-like,  medium  on  the  local 
quantum  system  of  interest. 

A  simple  correction  to  this  mean  field  correlation  function,  within  the 
spirit  of  QC  dynamics,  is  to  run  trajectories  for  the  classically  treated  degrees  of 
freedom  which  give  rise  to  a  time  dependent  potential  for  the  quantum  degrees  of 
freedom.  That  is,  He(r,p;R)  in  Eq.  (II.3)  becomes  time  dependent  due  to  the  time 
dependence  of  the  classical  position,  R(t),  in  the  potential  energy  operator 
component  of  Hg.  Only  the  quantum  component  of  the  Hartree  correlation 
function  of  Eq.(II.12)  is  kept,  giving 

Cqct^)  =  exp(iEgt/h)  <bp(0)  I  <|>p(t)) ,  (11.30) 

which  we  call  the  quantum-classical  trajectory  approximation  (QCT)  for  the 
spectral  correlation  function. 

How  to  choose  the  classical  trajectories  in  the  general  case  of  any  QC 
mixed  method  is  problematic.  By  comparing  to  the  tdH  theory,  Eq.(II.9-ll),  one 
finds  that  averaging  over  a  bundle  of  classical  trajectories  maintains  the  spirit  of 
integrating  over  all  R  in  the  tdH  equations.^  This  scheme  is  self-consistent  and 
energy  conserving  if  the  classical  trajectories  respond  to  a  mean  force  (averaged 
over  the  quantum  wave  function)  as  given  by  Eq.(II.16-19)  above.  For  the  set  of 
initial  conditions  specified  in  the  paragraphs  following  Eq.(II.19,\  the  bundle  of 
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trajectories  collapses  to  a  single  trajectory  that  is  uniquely  determined  by  the  set  of 
forces  impinging  on  the  frozen  cluster  configuration  due  to  the  excitation  of  the 
electron.  Because  of  this  unique  initial  condition  giving  rise  to  an  unambiguous 
mean  trajectory  approximation,  it  is  convenient  to  compare  the  results  of  the  QCT 
theory  to  the  QGWP  theory  for  the  spectroscopy  applications  examined  below. 

We  note  that  compared  to  the  tdH  correlation  function  in  Eq.(III.12),  the 
trajectory  approximation  correlation  function  in  Eq.(II.30)  lacks  the  factor  < 
Xg(0)  I  Xg(t) )  involving  wave  functions  for  the  classical-like  degrees  of  freedom 
(and  it  also  lacks  the  Hartree  phase  term).  The  underlying  assumption  in  the 
trajectory  approximation  is  that  the  time  dependent  potential,  arising  due  to  R(t) , 
and  appearing  in  the  quantum  Hamiltonian  will  act  to  remedy  this  deficiency.  In 
nearly  rigid  systems  where  the  static  mean  field  approximation  is  fairly  good,  this 
seems  reasonable.  However,  in  more  supple  systems  where  energy  transfer 
between  the  quantum  and  classical  degrees  of  freedom  distorts  the  classical 
positions  significantly,  it  is  obvious  from  the  form  of  the  tdH  correlation  function 
in  Eq.  (III.  12)  that  the  overlap  of  ( Xg(0)*  XgW  )  is  indispensable  in  obtaining  the 
correct  spectrum.  For  example,  if  the  system  is  floppy  enough  that  (  Xg(0)  IXglt) ) 
goes  to  zero  at  some  time  t,  then  the  overall  correlation  function  is  zero  and  in 
general  a  time  dependence  in  the  quantum  Hamiltonian  will  not  be  able  to 
account  for  this. 

ELF  Time-energy  resolution 

Heller  has  pointed  out^31^  the  useful  connection  between  the  resolution  of 
the  absorption  spectrum  and  the  length  of  time  in  which  the  promoted  state 
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explores  the  potential  energy  surface.  Engel,  Schinke,  Henning  and  Metiu^30^ 
have  discussed  in  detail  the  method  that  we  use  to  obtain  dynamic  information  by 
taking  the  high  resolution  spectrum  and  smearing  it  out  to  generate  a  sequence  of 
low  resolution  spectra. 

To  generate  this  spectral  resolution  we  replace  Eq.  (II.4)  by 


o(co;x)  ~  co  Re  [  dt  exp(icot)  fT(t)  C(t).  (11.31) 

JQ 

The  "window"  function 

£c(t)  =  exp[-(t/x)2].  (11.32) 

cuts  off  the  time  evolution  of  the  overlap  integral  for  times  substantially  longer 
than  x.  In  the  frequency  domain,  this  window  function  smears  out  the  spectrum 
and  generates  a  low  resolution  version  having  the  resolution  on  the  order 
Aco=27c/t.  Thus,  a  spectrum  taken  with  the  resolution  Aco  contains  information 
about  the  motion  of  the  promoted  state  for  a  time  x  =  2x/Aco. 

The  introduction  of  this  resolution  time  is,  of  course,  an  arbitrary 
parameter  from  a  dynamics  perspective.  One  advantage  of  including  a  "classical 
wave  function"  is  that  vibrationally  broadened  line  widths  of  electronic  transitions 
can  develop  more  naturally  on  a  time  scale  determined  by  the  dynamics  of  the 
problem.  As  the  classical  wave  functions  move  away  from  their  initial  positions, 
the  overlap  of  the  classical-like  term,  (xg(0)  I X^))  >  in  Eq.(II.12)  decays  with  time. 
If  the  configuration  of  the  promoted  system  is  such  as  to  prevent  further 
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recurrences  of  <xg(0)  I  xg(t))  i*1  time  then  a  broaden  line  will  remain  in  the 
spectrum  that  will  be  independent  of  x  (for  x  large).  If  the  configuration  of  the 
promoted  system  is  such  as  to  allow  for  recurrences  of  (Xg(0)  f  Xg^))  in  time,  then 
vibrational  fine  structure  will  appear  as  x  increases  and  the  spectral  lines  will 
continue  to  narrow  with  increases  in  the  resolution  time  x  . 
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IEL  TWO  COORDINATE  TEST  CASE 

We  examine  here  the  electric  dipole  absorption  spectrum  of  a  colinear 
alkali-helium  dimer.  Since  this  two-coordinate  model  problem  can  be  solved 
numerically  at  both  the  full  quantum  mechanical  level  and  the  full  quantum 
Hartree  level,  we  can  evaluate  the  accuracy  of  the  QGWP  method.  Computational 
details  for  the  two-coordinate  test  case  are  given  in  Appendix  A. 

HL1  The  model 

The  model  described  in  Ref.  21  is  used  here.  We  consider  an  electron 
(located  at  r),  an  ion  with  charge  +1  (of  infinite  mass  at  the  coordinate  origin),  and 
a  helium  mass  atom  (at  R  ).  The  model  consists  of  two  body  interactions  among 
the  three  species. 

For  the  alkali  ion  -  electron  interaction  energy  Vei,  we  use 
Vei  (r)  =  -e2  /rc  for  r  <  rc 

=  -e2/r  for  rc<r  (III.  1) 

In  Eq.(III.l)  e  is  the  electron  charge  and  rc  is  a  cutoff  parameter  that  will 
be  specified  below.  The  expression  in  Eq.(III.l)  has  been  used  to  represent  alkali 
atoms  by  Parinello  and  Rahman^34^  ,  Selloni,Carnevali,Carr  and  Parinello^10^  , 
and  Haug  and  Metiu^17^  in  a  variety  of  three-dimensional  calculations. 

The  electron- helium  interaction  energy  is  that  used  by  Coker,  Berne  and 
Thirumalai  ^  (for  electron-xenon)  with  a  large  distance  cut  off, 

Veh(X)  =  X  fa/X4)  [  b/(c+X6)  -  1  ] 


(III. 2) 
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where 

X=  min{lr-Ri  ,  Rc  },  (III. 3) 

a  =  26.86  eV  A4  ,  b  =  108.0  A6  ,  c  =  83.29A6.  The  cutoff  distance  Rc  =  0.5A  is 
introduced  to  make  the  potential  finite  (and  constant)  at  short  distances.  The  total 
electron  interaction  energy  of  Eq.(II.3)  is  the  sum  of  these  two  terms 

Vea(r,R)  =  Vei(r)  +'  Veh(r,R)  .  (III.4) 

The  alkali  ion  -  helium  interaction  energy,  constituting  Va(R)  of  Eq.(II.2), 
is  taken  to  be  harmonic  with  co  =  O.leV  and  which  roughly  approximates  the 
realistic  potassium  ion  -  xenon  interaction  given  by  Budenholzer,  Gislason  and 
Jorgensen^. 

ITI.2  Absorption  spectrum 

In  Figure  1,  we  present  the  fundamental  transition  of  the  electronic 
absorption  spectrum  for  this  model  problem  calculated  with  the  quantum  Hartree 
approximation,  Fig. la,  and  the  QGWP  method,  Fig.lb.  The  quantum  Hartree 
method  is  calculated  using  Eqs.(II.9-12)  with  quantum  dynamics  for  both  degrees 
of  freedom.  This  Hartree  method  produces  a  spectrum  in  which  the  major  peaks 
are  essentially  exact  when  compared  to  a  full  quantum  calculation^21^.  The 
fundamental  electronic  transition  given  in  Fig.  1  shows  the  vibrational  fine 
structure  arising  from  the  O.leV  dimer  vibration  in  the  3  split  lines.  Ref.  21 
presents  details  regarding  the  small  spurious  bumps  and  dips  which  appear 
surrounding  the  major  peaks  which  are  related  to  the  lack  of  quantum 
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correlations  inherent  in  a  Hartree  product  assumption.  Since  the  QGWP 
approximation  uses  a  Hartree  product,  we  cannot  expect  it  to  do  better  than  time 
dependent  Hartree. 

As  seen  in  Fig.  1,  the  QGWP  spectrum  is  a  reasonable  representation  of 
the  quantum  Hartree  spectrum.  The  fundamental  electronic  transition  appears 
with  vibrational  fine  structure  at  nearly  the  same  frequencies  and  with  relative 
amplitudes  preserved  as  well.  The  width  of  the  spectral  lines  in  Fig.l  is 
determined  by  the  cutoff  time  of  Eq.(II.32)  ;  a  choice  of  x=100fs  was  used  in  this 
case.  This  result  provides  us  with  a  measure  of  confidence  in  applying  the  QGWP 
method  to  more  complex  and  realistic  situations  where  a  comparison  to  a  more 
precise  method  is  computationally  out  of  hand. 
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IV.  KXen  CLUSTER 

We  now  apply  and  compare  the  QCT  and  the  QGWP  methods  to  a  more 
complex  and  physically  realistic  model.  We  want  to  examine  how  the  the 
absorption  spectrum  of  the  K  atom  is  modified  by  the  interaction  with  the  Xe 
atoms  forming  the  cluster.  This  difficult  physical  problem  can  be  examined 
satisfactorily  if  we  work  with  an  atom  consisting  of  one  valence  electron 
interacting  with  a  simple  two  body  pseudopotential  centered  on  the  ion.  The 
potential  is  spherically  symmetric,  energy  independent,  and  local.  It  is  further 
assumed  that  electron-solvent  and  the  ion-solvent  interactions  are  two  body 
interactions.  While  some  of  these  assumptions  can  be  improved,  we  see  no  reason 
at  this  time  to  increase  the  complexity  of  the  calculations  by  treating  some  of  the 
interactions  at  a  higher  level  of  accuracy  than  others.  Our  purpose  here  is  to 
examine  the  methodology.  In  Section  IV.  1  we  give  a  brief  summary  of  the 
potentials  used.  The  computational  details  for  this  application  are  given  in 
Appendix  A. 

IV.  1  Interaction  Potentials 

For  the  ion  -  electron  interaction  energy,  Vei ,  we  use  Eq.(III.l)  with 
rc  =2.2 1A  which  results  in  the  correct  ionization  potential  for  the  K  atom  (4.34eV). 
We  note  that  the  spectrum  will  be  distorted  in  the  high  frequency  region 
corresponding  to  transitions  to  final  states  whose  outer  turning  points  are  close  to 
or  exceed  rc.  In  spite  of  the  simplicity  of  this  psuedopotential,  the  low  resolution 
absorption  spectrum  of  an  isolated  K  atom,  calculated  with  this  potential,  is 
similar  to  the  measured  one  (see  Table  I,  to  be  discussed  below). 

The  xenon  -  electron  interaction  energy  is  that  used  by  Coker,  Berne  and 
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Thirumalai  with  a  large  distance  cutoff.  The  total  electron  interaction  with 
the  N  xenon  (solvent)  atoms  is 

Ves  =  X  (a^n4)  C  b/(c+Rn6)  - 1  ]  (IV. 2) 

n 

where 

Rn  =  min  {  X„  ,  Rc  },  (IV.3) 

is  the  distance  between  the  electron  and  the  n-th  atom,  a  =  26.86  eV  A4  ,  b  = 
108.0  A6  ,  and  c  =  83.29A6.  The  cutoff  distance  Rc  is  introduced  to  make  the 
potential  finite  (and  constant)  at  short  distances.  The  magnitude  of  Rc  determines 
the  "transparency"  of  the  xenon  atoms  to  the  electron  and  we  set  Rc  =  0.5A.  The 
effect  of  varying  Rc  is  examined  in  Ref.  17.  This  potential  has  a  shallow  attractive 
well  of  0.48eV  at  2.4  A,  which  for  a  single  Xe  cannot  bind  the  electron. 

The  total  electron-atom  potential  used  in  Eq.(II.3)  is  the  sum  of  these 

terms 

Vea(r,R)  =  Vei  +  Ves .  (IV. 4) 

For  the  xenon-xenon  interaction,  we  use  a  pairwise  Lennard-Jones 
potential^, 

vss  =  (1/2)2  4e[  •  WTW6  +  (°/R„m)12  3.  (IV. 5) 

n,m 

with  the  parameters  e  =  1.9733  x  10'2eV  and  a  -  4.055 lA.  Rnm  is  the  distance 
between  the  n-th  and  m-th  xenon  atoms. 

The  potassium  ion  -xenon  interaction  energy  is^36^ 
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Vis  =  I  Cg/R^3  -  C6  /  Rjn6  -  C4  /  R^4 ,  (IV.6) 

n 

where  Rin  is  the  distance  between  the  n-th  xenon  atom  and  the  potassium  ion, 

C9=  9.70  x  10+3  eV  A9  ,C6=  1.4104  x  10+2eV  A6  and  C4=  2.908  x  10+1  eV  A4  .  This 

* 

potential  has  a  weak  attractive  well  of0.15eV  at  a  distance  of  3.4A. 

The  total  inter-atomic  interaction  used  in  Eq.(II.2)  is  then 

va(R)  =  Vis  +  vss . 


(IV. 7) 
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IV.2  The  spectrum  of  the  isolated  potassium  atom 

The  absorption  spectrum  of  the  isolated  K  atom  is  obtained  by  including  in 
the  Hamiltonian  only  the  ion  -  electron  interaction  in  Eq.(IV.4).  The  absorption 
spectrum  is  given  in  Fig.  2a.  We  have  examined  it  to  provide  a  reference  for 
determining  the  spectral  changes  caused  by  the  presence  of  the  Xe  atoms.  The 
calculated  energies,  oscillator  strengths  and  peak  frequencies  for  the  K  atom 
using  the  computational  parameters  discussed  in  Section  III  are  compared  to 
measured  values^37^  in  Table  I.  The  frequency  of  the  fundamental  peak  and  the 
oscillator  strength  of  the  first  three  peaks  in  the  calculated  spectrum  are  close  to 
the  experimental  values  (see  rows  1-2  of  Table  I  for  the  oscillator  strengths  and 

the  peak  frequency  in  our  model  of  K  versus  real  K). 

The  first  two  peaks  in  the  spectrum  are  rather  similar  to  those  measured 
for  the  K  atom  and  our  results  concerning  these  peaks  are  relevant  to 
experiments  performed  on  the  KXen  system.  The  calculated  electronic 
transitions  are  found  at  about  1.64,  3.08  and  3.49eV  while  experimentally^3^  the 
transitions  occur  at  1.61,  2.67  and  3.80  eV.  Errors  in  the  higher  energy  peaks  are 
expected  since  the  absorption  energies  of  these  peaks  are  in  the  energy  range 
where  the  computational  grid  is  cut  off. 

IV. 3  KXes  cluster 

The  KXeg  cluster,  for  the  potentials  specified  in  Section  III,  has  a  ground 
state  minimum  energy  configuration  of  octahedral  symmetry  with  K  at  the  center 
and  the  Xe  atoms  located  3.45A  from  the  K  center.  The  electronic  absorption 
spectrum  of  this  cluster  is  calculated  using  three  alternative  dynamics 
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approximations.  In  the  first  case,  the  cluster  is  frozen  in  its  ground  state 
equilibrium  configuration  and  the  spectrum  is  calculated  by  letting  the  promoted 
wave  function  interact  with  the  Xe  atoms.  In  the  second  and  third  cases  we  allow 
the  atoms  to  move  in  response  to  the  excitation  of  the  electron,  where  this  motion 
is  calculated  by  the  QCT  method  and  also  by  the  QGWP  method. 

We  place  the  six  Xe  atoms  symmetrically  along  the  x,y  and  z  axis  with  the 
K  atom  at  the  origin.  The  spectra  shown  in  Fig.  2  were  calculated  with  an  electric 
field  polarization  direction  in  the  {  0,0,1  )  direction  (along  one  of  the  K-Xe  axes)  ; 
and  the  spectrum  in  Fig.  3  uses  an  electric  field  polarization  direction  in  the  ( 
1,1,1  )  direction  (between  the  K-Xe  axes). 

In  Fig.  2b  we  present  the  spectrum  using  the  frozen  KXeg  cluster.  The 
spectral  lines  are  shifted  due  to  the  interaction  of  the  electron  with  the 
surrounding  Xe  atoms.  A  resolution  time  of  t=150?s  is  used  for  the  spectre1 
window  function,  in  Eq.(II.32).  This  frozen  cluster  spectrum  has  a  blue  shift  of 
about  0.33eV  compared  to  the  lone  K  atom  spectrum.  The  electronic  transitions 
are  at  about  1.97,  3.36  and  4.13eV.  The  relative  intensities  of  the  first  few  lines  are 
changed  only  slightly  by  the  presence  of  the  Xe  (see  Table  I,  row  2  versus  row  3). 

In  Fig.  2c  we  present  the  spectrum  of  the  KXeg  cluster  using  the  more 
traditional  quantum-classical  trajectory  method  (QCT  of  Section  II. E)  in  which 
the  nuclear  overlap  is  left  out.  A  resolution  time  of  t=150fs  is  used  for  the  spectral 
window  function,  Eq.(II.32)  ,  for  this  case.  The  spectrum  for  the  QCT  method 
shows  a  dramatic  anomaly  on  the  red  side  of  the  fundamental  transition  which 
arises  as  the  cluster  relaxes,  in  response  to  the  promotion  of  the  electron.  This 
anomaly  begins  to  appear  at  a  resolution  time  of  about  t  =  60fs  and  becomes  more 
pronounced  with  increas. ng  values  of  t  .  For  values  of  t  >  200fs,  the  anomaly 
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dominates  the  spectrum.  We  point  out  that  this  ringing  anomaly  is  not  an  artifact 
caused  by  aliasing  in  the  Fourier  transforms,  but  is  a  real  effect  due  to  the 
incorrect  feedback  between  the  quantum  and  classical  degrees  of  freedom  The 
relative  peak  positions  (  at  about  1.96,  3.36  and  4.13eV  )  and  intensities  of  the  lines 
in  the  QCT  spectrum  in  Fig.  2c  are  both  in  agreement  with  the  frozen  cluster 
spectrum  of  Fig.  2b  (see  also  Table  I,  rows  3  and  4). 

In  Fig.  2d,  the  spectrum  of  the  KXe6  cluster  using  the  QGWP  method  is 
presented.  In  this  cluster  the  motion  of  the  classical  wave  function,  in  the  QGWP 
approximation,  provides  an  intrinsic  decay  for  the  correlation  function  and  we  do 
not  ne^d  to  use  a  window  function.  The  line  widths  are  "natural"  in  the  sense 
that  they  represent  the  decay  of  the  vibrational  overlap,  which  takes  place  in  about 
150fs. 

The  GWP  wave  function  ,in  Eq.(II.25),  separates  various  physical  effects 
entering  in  the  nuclear  wave  function.  It  contains  an  amplitude  (  which  for 
simplicity  we  keep  constant  in  width)  ,and  a  phase  given  by  the  classical  action 
along  the  trajectory  traced  by  the  center  of  the  packet  along  with  a  factor  exp( 
ipt»(r  -  rt)/  h  ]  associated  with  the  momentum  distribution  of  the  packet.  The 
spectrum  is  practically  unchanged  if  we  remove  the  phase  of  the  GWP  wave 
function.  The  main  contribution  to  the  spectrum  comes  from  the  amplitude  term, 
which  causes  a  decay  of  the  nuclear  overlap  entering  in  the  cross  section  formula. 
In  fact  if  we  suppress  the  phase  information  of  the  classical  wave  function  for  this 
system  and  calculate  a  spectrum  for  comparison  with  the  spectrum  of  Fig.  2d  ,  we 
find  marginally  perceptible  differences  on  the  scale  give  in  Fig.  2d. 

From  the  comparison  of  Fig.  2c  and  2d,  it  is  clear  that  the  QCT 


28 


approximation  has  greater  trouble  with  negative  cross  section  values  than 
QGWP.  Never-the-less  small  negative  cross  sections  do  occur  in  the  QGWP 
spectrum.  Since  such  small  anomalies  also  occur  in  the  exact  Hartree 
calculations  (  see  Fig.  1  above  and  Ref.  16  and  21),  it  is  quite  possible  that  the 
anomalies  seen  in  Fig.  2d  are  due  to  the  Hartree  product  assumption  made  in 
Eq.(II.8)  and  that  no  method  using  a  Hartree  product  wave  function  would 
improve  on  this. 

Two  other  points  need  to  be  made  when  the  QCT  and  the  QGWP  methods 
are  compared.  The  relative  oscillator  strength  of  the  spectral  lines  given  by  the 
two  methods  are  considerably  different  (  see  Table  I).  The  relative  strength  of  the 
2p,3p,4p,...  transitions  in  the  series,  is  principally  governed  by  the  electronic 
time-correlation  matrix  element,  <<{>p(0)  1  <bp(t»  in  Eq.(II.12),  which  does  not  vary 
between  the  two  methods.  However  the  Hartree  phase  term,  involving  q(t)  in 
Eq.(II.12),  which  is  left  out  in  the  QCT  method  causes  most  of  the  intensity 
difference  between  the  two  methods.  This  can  be  seen  from  Table  I  by  comparing 
line  5  (the  QCT  method)  ,  with  line  6  (  QCT  with  the  Hartree  phase  term  added  ), 
and  line  7  (the  QGWP  method) . 

A  second  difference  between  the  QCT  and  the  QGWP  methods  is  the 
energy  shift  in  the  spectrum.  The  QGWP  electronic  transitions  are  at  about  1.39, 
2.79  and  3.55eV  which  is  a  red  shift  of  about  0.57eV  (  for  the  fundamental)  with 
respect  to  the  QCT  method.  This  shift  is  readily  traced  to  the  phase  information  in 
the  classical  wave  function  and  the  Hartree  phase  terms.  The  QGWP  spectrum  is 
also  red  shifted  by  about  0.25eV  compared  to  the  spectrum  of  the  isolated  K  atom. 
While  we  are  not  aware  of  spectral  data  for  small  K/Xe  clusters,  experimental 
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data  for  K  defects  in  Xe  matrices  at  1.54°K  ^38c^  suggests  three  different  sites  for  K 
in  Xe,  each  of  which  show  triplet  patterns.  Two  sites  have  centers  (of  the  triplet) 
showing  blue  shifts  of  the  fundamental  K  absorption  line  of  about  0.05eV  and 
0.02eV  while  the  third  shows  a  red  shift  of  about  O.OleV.  Assuming  that  our 
simplified  potentials  are  accurate,  the  shift  difference  that  we  find  could  be 
explained  by  noting  that  the  Xe  matrices  do  not  have  the  flexibility  of  the  KXe6 
cluster  where  the  Xe  atoms  have  considerable  room  to  relax.  Over  the  course  of 
the  150fs  decay  time  of  the  nuclear  overlap  of  the  Xe  atoms, the  Xe  aligned  with  the 
electric  field  polarization  (  which  are  most  strongly  affected  by  the  promoted 
electron)  relax  by  about  0.085  Angstroms,  while  the  Xe  atoms  perpendicular  to  the 
electric  field  direction  relax  by  about  0.039  Angstroms.  Such  a  large  relaxation  is 
less  likely  in  a  matrix  and  this  difference  may  explain  the  difference  in  the 
spectral  shift.  The  Xe  atoms  in  this  small  cluster  tend  to  accommodate  the 
"needs"  of  the  electron  and  are  therefore  displaced  in  a  way  that  lowers  the  energy 
of  the  excited  state. 

As  noted  above,  in  the  K/Xe  matrix  experiments^38*3’0^  the  fundamental  K 
line  splits  into  a  triplet  separated  by  about  (0.01-0.02)eV.  The  first  excited  state  of 
the  valence  electron  in  the  octahedral  KXe6  cluster  is  also  three-fold  degenerate. 
Due  to  the  symmetry  of  the  cluster  and  the  polarization  directions  chosen,  this 
triplet  does  not  appear  in  Figs.  2  and  3.  For  example,  in  the  case  of  the  (  0,0.1  > 
polarization  direction,  although  the  octahedral  symmetry  is  broken  by  a 
tetragonal  elongation  (  ->  D4h  )^39^  as  the  Xe  atoms  along  the  polarization 

direction  back  away  from  the  K  center  when  the  electron  is  promoted,  the  2-fold 
states  of  the  doublet  are  still  orthogonal  to  the  polarization  direction.  In  the  case 
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of  the  <  1,1,1  )  polarization  direction,  the  degeneracy  is  not  lifted.  For  arbitrary, 
polarization  directions  the  splitting  may  appear  depending  on  the  magnitude  and 
time  scale  of  the  configurational  changes  occurring  in  the  cluster  upon  excitation 
(see  Section  II.F  above). 

The  normal  mode  frequencies  of  the  excited  cluster  range  from  1.4  x  10'3 
to  2.5  x  10‘2  fs'1,  giving  recurrence  times  in  the  range  of  250  *  4500fs.  While  our 
main  interest  has  been  in  short  time  calculations,  we  have  also  carried  out  a 
calculation  meant  to  observe  the  appearance  of  vibrational  structure  in  the 
electronic  spectrum.  In  Fig.  3a  ,  the  spectrum  of  the  KXe6  cluster  using  the 
QGWP  method  is  presented  for  a  polarization  direction  chosen  in  the  (  1,1,1  ) 
direction  which  perturbs  all  six  Xe  atoms  equally.  For  this  polarization  direction 
the  anomalous  features  that  appeared  in  the  QCT  spectrum  in  Fig.  2c  are  again 
minimal.  Also  appearing  in  this  spectrum  is  the  beginning  cf  a  vibrational 
structure,  which  is  more  apparent  when  attention  is  focused  upon  Fig.  3b  in 
which  the  fundamental  peak  of  Fig.  3a  is  blown-up.  The  vibrational  structure  is 
characterized  by  an  energy  spacing  of  about  0.00166eV  (  13.4  cm'1 )  which  implies 
a  corresponding  time  periodicity  of  about  2500fs.  All  six  Xe  atoms  are 
dynamically  equivalent  for  this  polarization  direction  and  a  representative  Xe-K 
distance  is  displayed  in  Fig.  4a  for  4000fs.  The  motion  of  the  Xe  atoms  is  driven  by 
the  energy  flow  from  the  quantum  mechanical  degrees  of  freedom  (the  excited 
electron)  into  the  classical-like  degrees  of  freedom.  This  energy  flow  is  indicated 


in  Fig.4b  for  4000fs.  The  classical  energy  ,  denoted  ”C"  in  Fig.  4b,  is 

EC=I  P„2/2Mn  +  Va<R) 


(IV. 8) 


and  is  measured  along  the  left  axis  while  the  quantum  energy,  denoted  "Q"  in 
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Fig.  4b,  is 

Eq  =  <♦„  I  p2/2m  +  Vea(r.(R»  I  <(.p>  /  I  *,„>  (IV.9) 

and  is  measured  along  the  right  axis  of  the  figure.  From  both  Fig.  4a  and  4b  a 
vibrational  time  cycle  of  about  2500fs  is  apparent.  Note  also  that  the  vibrational 
mode  is  not  just  excited  in  this  case,  but  rather  that  energy  flows  into  the  classical 
coordinates  over  the  first  700fs  and  then  drains  out  of  these  coordinates  in  the 
1700-2400fs  time  span,  after  which  it  cycles  back  in  again.  This  emphasizes  the 
significant  coupling  of  the  electronic  and  vibrational  degrees  of  freedom  in  this 
system. 

We  note  that  for  the  calculations  in  Fig.  4  the  total  energy  of  the  system, 


quantum  plus  classical,  is  conserved  to  better  than  0.7  percent  over  the  4000fs  run 
time.  Further  computational  details  are  given  in  Appendix  A. 
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V.  Summary 

We  have  examined  a  computational  method  in  which  a  set  of  classical- 
like  degrees  of  freedom  are  described  by  Gaussian  wave  packets  while  the  others 
are  treated  quantum  mechanically.  This  method  was  applied  to  a  two-coordinate 
model  problem  in  which  exact  numerical  results  could  be  obtained  and  we  found 
that  the  method  works  fairly  well.  We  have  also  calculated  the  electronic 
absorption  spectrum  for  a  KXeg  cluster  in  a  frequency  range  where  the  valence 
electron  of  the  K  atom  is  excited.  The  method  does  not  suffer  from  the  dramatic 
failure  (  negative  absorption  cross  section)  seen  when  Xe  motion  is  treated 
classically.  Moreover,  it  generates  a  reasonable  spectrum.  The  use  of  Gaussian 
wave  packets  (GWP)  is  justified  by  the  observation  that  the  overlap  of  the  nuclear 
wave  function  is  likely  to  decay  very  rapidly,  requiring  only  a  short  time 
calculation,  for  which  the  GWP  method  is  accurate.  The  decay  time  is  shorter  for 
systems  with  many  degrees  of  freedom  and  the  accuracy  of  the  method  increases 
as  the  system  becomes  more  complete. 

The  most  interesting  applications  that  we  foresee  are  for  several  quantum 
systems  imbedded  in  a  condensed  medium  or  a  large  cluster.  In  such  systems 
the  width  of  the  spectrum  is  due  to  the  vibrations  of  the  medium  as  well  as 
inhomogeneous  broadening.  This  method  can  calculate  the  vibrational  width  and 
the  vibrational  structure. 
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Appendix  A:  Computational  details 

The  initial  ground  state  wave  functions  for  both  the  quantum  and  Hartree 

cases  are  calculated  by  acting  with  exp(-aH)  on  an  arbitrary  trial  wave  function 

with  the  appropriate  symmetry.  As  the  real  number  a  becomes  larger,  the  wave 

function  converges  towards  the  ground  state  of  the  Hamiltonian  H.  For 

propagation,  we  use  a  method  proposed  by  Fleck,  Morris  and  FeiJ33l  The 

numerical  procedure  has  been  described  by  Hellsing,  Nitzan  and  Metiu^40^  and 

has  been  used  by  others  in  similar  contexts^4 

The  quantum  time  propagation  is  performed  with  the  algorithm  of  Fleck, 

Morris  and  Feit^33^,  as  discussed  in  Ref.  17.  In  the  QC  calculations,  we  propagate 

the  electron  wave  function  for  one  time  step,  calculate  the  mean  potential,  update 

the  classical  positions  ,  update  the  classical  coordinates  in  the  quantum  electron 

Hamiltonian,  propagate  the  quantum  wave  function  for  another  d me  step,  etc.-  '- 

In  the  two-coordinate  model  for  the  quantum  calculations  ,  a  time  step  of 

O.Olfs  is  used  and  a  spatial  grid  having  64  points  and  a  spacing  of  0.6  A  in  the  r- 

coordinate  and  a  spatial  grid  having  96  points  with  a  spacing  of  0.067 A  in  the  R- 

coordinate.  The  results  were  checked  for  convergence  with  respect  to  grid 

spacing  and  grid  size.  The  classical  trajectories  are  generated  by  using  the  Verlet 

algorithm^42^  with  a  time  step  of  0.01  fs. 

In  the  KXeg  model- for  the  quantum  coordinates,  a  time  step  of  O.Olfs  is 

« 

used  and  a  3-D  spatial  grid  having  36  points  and  a  spacing  of  0.7  A  in  each 
direction.  The  classical  trajectories  are  generated  by  using  the  Verlet  algorithm 
with  a  time  step  of  0.05  fs.  The  results  were  checked  for  convergence  with  respect 
to  grid  spacing,  grid  size,  and  time  step. 


34 


ACKNOWLEDGMENT 

This  work  is  supported  by  NSF  CHE91-12926.  Acknowledgment  is  also 
made  to  the  donors  of  The  Petroleum  Research  Fund,  administered  by  the  ACS, 
for  partial  support  of  this  research  and  also  to  a  grant  from  Lehigh  University. 


35 


REFERENCES 

[1]  The  number  of  articles  applying  such  a  technique  is  so  large  that  we 

mention  here  only  a  few  recent  articles  [2-21]  which  illustrate  the  procedure 

and/or  which  are  specifically  discussed  in  this  article. 

[2]  J.T.  Muckerman,  I.  Rusinek,  R.E.  Roberts  and  M.  Alexander,  J.  Chem. 
Phys.  65,  2416(1976). 

[3]  W.H.  Miller,  J.  Chem.  Phys.  68,  4431  (1978) . 

[4]  S.D.  Augustin  and  H.  Rabitz,  J.  Chem.  Phys.  69,  4195  (1978). 

[5]  D.J.  Diestler,  J.  Chem.  Phys.  78,  2240  (1983). 

[6]  J.  Olsen  and  D.  Micha,  Int.  J.  Quantum  Chem.  22,  971  (1982). 

[7]  G.D.  Billing,  Chem.  Phys.  70,  223  (1982). 

[8]  A.M.  Roberts  and  A.E.  DePristo,  Surf.  Sci.  134,  338  (1983). 

[9]  Carr  and  Parrinello,  Phys.  Rev.  Lett.  55,  2473  (1985). 

[10]  A.  Selloni,  P.  Carnevali,  R.  Car  and  M.  Parrinello,  Phys.  Rev.  Lett.  59, 
823  (1987). 

[11]  G.  Wahnstrom,  B.  Carmeli  and  H.  Metiu,  J.  Chem.  Phys.  88,  2478  (1988) 

[12]  C.Z.  Wang,  C.T.  Chan  and  KM.  Ho,  Phys.  Rev.  B  39,  8586  (1989). 

[13]  R.  Alimi,  R.  B.  Gerber,  A.D.  Hammerich,  R.  KoslofF  and  M.  A.  Ratner  ,  J. 
Chem.  Phys.  93,  6484  (1990). 

[14]  F.  Webster,  P.J.  Rossky,  and  R.A.  Friesner,  Comp.  Phys.  Comm.  63,  493 
(1991). 

[15]  R.N.  Barnett,  U.  Landman  and  A.  Nitzan,  J.  Chem  Phys.  94,  608  (1991). 

[16]  E.J.  Heller,  J.  Chem.  Phys.  94,  2723  (1991). 

[17]  K.  Haug  and  H.  Metiu,  J.  Chem.  Phys.  95,  5670  (1991). 

[18]  B.  Hartke  and  E.A.  Carter,  Chem.  Phys.  Lett.  189,  358  (1992) 

[19]  M.J.  Field,  J.  Chem.  Phys.  96,  4583  (1992). 

[20]  E.  Deumenes,  A.  Diz,  H.  Taylor  and  Y.  Ohrn,  J.  Chem.  Phys.  96,  6820 
(1992). 


36 


[21]  K.  Haug  and  H.  Metiu,  J.  Chem.  Phys.  97,  4781  (1992). 

[22]  J.B.  Delos,  W.R.  Thorson  and  S.K.  Knudson,  Phys.  Rev.  A  6,  709  (1972);  J.B. 
Delos  and  W.R.  Thorson  Phys.  Rev.  A  6,  720  (1972). 

[23]  H.D.  Meyer  and  W.  H.  Miller,  J.  Chem.  Phys.  70,  3214  (1979). 

[24]  R.  B.  Gerber,  V.  Buch  and  M.  A.  Ratner,  J.  Chem.  Phys.  77,  3022  (1982). 

[25]  D.  Thirumalai,  A.J.  Bruskin  and  B.J.  Berne,  J.  Chem.  Phys.  83,  230  (1985). 

[26]  N.  Makri  and  W.H.  Miller,  J.  Chem.  Phys.  87,  5781  (1987). 

[27]  Z.  Kotler,  A.  Nitzan,  and  R.  Kosioff,  Chem.  Phys.  Lett,  153,  483  (1988). 

[28]  A.  Garcia-Vela  and  R.  B.  Gerber,  J.  Chem.  Phys.  98,  427  (1993). 

[29]  (a)  E.J.  Heller,  J.  Chem.  Phys.  62,  1544  (1975);  65,  1289  (1976);  65,  4979 
(1976);  67,3339(1977). 

(b)  S.  Sawada,  R.  Heather,  B.  Jackson  and  H.  Metiu,  J.  Chem.  Phys.  83, 
3009  (1985);  R.  Heather  and  H.  Metiu,  J.  Chem.  Phys.  84,  3250  (1986). 

[30]  V.  Engel,  R.  Schinke,  S.  Hennig  and  H.  Metiu,  J.  Chem.  Phys.  92,  1  (1990). 

[31]  E.J.  Heller,  J.  Chem.  Phys.  68,  2066  (1978);  68,  3891  (1978). 

[32]  (a)  P.A.M.  Dirac,  Proc.  Camb.  Phil.  Soc.  25,  62  (1929). 

(b)  A.D.  McLachlan,  Mol.  Phys.  8,  39  (1964). 

[33]  J.A.  Fleck,  Jr.,  J.R.  Morris  and  M.D.  Feit,  Appl.  Phys.  10,  129  (1976). 

[34]  M.  Parrinello  and  A.  Rahman,  J.  Chem.  Phys.  80,860  (1984). 

[35]  D.F.  Coker,  B.J.  Berne  and  D.  Thirumalai  ,  J.  Chem.  Phys  86,  5689  (1987). 

[36]  F.E.  Budenholzer,  E.A.  Gislason,  A.D.  Jorgensen,  J.  Chem.  Phys.  66, 

4832  (1977). 

[37]  (a)  H.G.  Kuhn,  Atomic  Spectra  (Longmans,  Green  and  Co.  ,  1969). 

(b)  A.R.  Striganov  and  N.S.  Sventitskii,  Tables  of  Spectral  Lines  of 
Neutral  and  Ionized  Atoms,  (  Plenum  Press,  New  York,  1968). 


37 


[38]  For  an  overall  view,  see:  (a)  D.W.  Ball,  Z.H.  Kafafi,  L.  Fredin,  R.H.  Hauge 
and  J.L.  Margrave,  editors,  A  Bibliography  of  Matrix  Isolation 
Spectroscopy  1954-1985  (Rice  University  Press,  Houston,  1988). 

The  spectroscopy  of  alkali  atoms  in  rare  gas  matrices 
can  be  found  going  back  to  : 

(b)  M.  McCarty,  Jr.  and  G.W.  Robinson,  Mol.  Phys.  2,  415  (1959);  W. 
Weyhmann  and  F.M.  Pipkin,  Phys.  Rev.  137,  A490  (1965);  B.  Meyer,  J. 
Chem.  Phys.  43,  2986  (1965); 

and  a  sample  of  more  recent  work  includes  : 

(c)  C.  Samat,  J.L.  Rose,  B.E.  Williamson,  and  P.N.  Schatz,  Chem.  Phys. 
Lett.  142, 557  (1987) ; 

(d)  M.E.  Fajardo,  J.  Chem.  Phys.  98,  110  (1993). 

[39]  F.A.  Cotton  and  G.  Wilkinson,  Advanced  Inorganic  Chemistry  ,  4th  ed.  , 
p.  51,  (John  Wiley  and  Sons,  New  York,  1980). 

[40]  B.  Hellsing,  A.  Nitzan  and  H.  Metiu  ,  Chem.  Phys.  Lett.  123,  523  (1986). 

[41]  J.  Schnitker,  K.A.  Motakabbir,  P.J.  Rossky  and  R.  Friesner,  Phys.  Rev  Lett. 
60,  456  (1988);  K.A.  Motakabbir  and  P.J.  Rossky,  Chem.  Phys.  129,  253 
(1989);  K.A.  Motakabbir,  J.  Schnitker  and  P.J.  Rossky,  J.  Chem.  Phys.  90, 
6916(1989). 

L.  Verlet,  Phys.  Rev.  159,  98  (1967). 


[42] 


38 


TABLE  I: 

The  oscillator  strength  (Fn)  for  the  first  three  peaks  in  the  spectrum  is  given  in 
columns  2-4.  The  ionization  potential  (IP  in  eV)  is  given  in  column  5.  The 
fundamental  peak  frequency  (co  ,  in  eV)  in  the  spectrum  is  listed  in  column  6. 


System 

n=l 

Fn  in  series 

n=2 

n=3 

IP 

(ev) 

peak 

co(eV) 

K  atom  (expt.)  a 

0.968 

0.009 

0.001 

4.34 

1.61 

K  atom  (calculated) 

0.988 

0.009 

0.001 

4.34 

1.64 

KXe6  (calculated): 

Frozen  Xe 

0.962 

0.027 

0.008 

4.77 

1.97 

QCT  method 

0.953 

0.027 

0.008 

4.77 

1.96 

QCT  method  +  r|(t) 

0.730 

0.159 

0.031 

4.77 

b 

QGWP  method 

0.755 

0.141 

0.027 

4.77 

1.39 

(a)  Best  fit  to  a  variety  of  experimental  and  theoretical  data.'37a^ 

(b)  The  peak  frequency  value  for  this  case  (7.94eV)  is  meaningless  since  energy 
shifts  are  being  included  via  the  Hartree  phase  term  which  are  not  being 
corrected  for  from  the  phase  of  the  classical  wave  function. 
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FIGURE  CAP!  IONS 

Fig.  1  The  absorption  spectrum  (in  arbitrary  units)  versus  frequency  (in  eV)  for 
the  two-coordinate  model  problem  using:  (a)  the  time-dependent  Hartree  method  ; 
(b)  the  QGWP  approximation. 

Fig.  2  The  absorption  spectrum  (in  arbitrary  units)  versus  frequency  (in  eV)  with 
the  (  0,0,1 )  polarization  direction  (  along  the  K-Xe  axis  )  for: 

(a)  the  lone  K  atom  ;  (b)  the  frozen  KXeg  cluster  ;  (c)  the  KXeg  cluster  using  the 
QCT  approximation  ;  (d)  the  KXeg  cluster  using  the  QGWP  approximation. 

Fig.  3  The  absorption  spectrum  (in  arbitrary  units)  versus  frequency  (in  eV)  with 
the  (  1,1,1  )  polarization  direction  (between  the  K-Xe  axes)  for  the  KXeg  cluster 
using  the  QGWP  approximation,  (a)  The  full  spectrum  on  the  0-4  eV  energy 
scale  ;  (b)  the  fundamental  transition  from  part  (a)  showing  the  vibrational  fine 
structure. 

Fig.  4  For  the  KXeg  cluster  described  in  Fig.  3:  (a)  a  representative  coordinate 
R(t)  ,  in  Angstroms,  giving  the  Xe-K  distance  for  4000fs  ;  (b)  curve  C  is  the  energy 
of  the  classical  degrees  of  freedom  ('n  eV  ,  on  the  left  axis)  and  curve  Q  is  the 
energy  of  the  quantum  degrees  of  freedom  (in  eV  ,  on  the  right  axis) 
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